function f = jacobNL(thetaNL,delta,Data)

% jacobNL - This function is calculate the derivative of the mean value
% delta wrt the parameter using the implicit function theorem (see Nevo
% http://faculty.wcas.northwestern.edu/~ane686/supplements/Ras_guide_appendix.pdf
% for the Nested Logit
% 
% jacobNL(thetaNL,delta,Data)
%
% Inputs:
%    input1 - Non linear parameters nested logit
%    input2 - Data
%
% Outputs:
%    Derivative of Delta wrt thetaNL
%
% Subfunctions: NLShareCalculation; multinv; multiprod
%

% Author: Laura Grigolon and Frank Verboven
% August 2012;

% Unpack
prods           = Data.prods;
nmkt            = Data.nmkt;
Nnodes          = Data.Nnodes;
qweightrprods   = Data.qweightrprods;
multidumseg     = Data.multidumseg;

xvuSEG          = Data.xvuSEG;
multixvuSEG     = Data.multixvuSEG;

xvuRC           = Data.xvuRC;
multixvuRC      = Data.multixvuRC;

% Market Shares
[sij,sijg,~,~,~, numer1,denom1,~,denom2] = NLShareCalculation(thetaNL,delta,Data);

if numel(thetaNL) == 3

    RCSeg           = thetaNL(1);
    RCCont          = thetaNL(2);
    Sigseg          = thetaNL(3);
else
    
    RCSeg           = 0; % this is for the correctly specified DGP
    RCCont          = thetaNL(1);
    Sigseg          = thetaNL(2);
end

% Derivative wrt delta

part1          = sum((1./(1-Sigseg) * sij),2);
part1          = bsxfun(@times,part1,eye(prods));          % diagonal in multiple dimensions

sijtransp      = permute(sij,[2 1 3 4]);
part2          = (Sigseg/(1-Sigseg))*multiprod(sijg,sijtransp);

part3          = multiprod(sum(sij,2),sum(sijtransp,1));

derShareDeltij = (part1 - part2 - part3);
derShareDelta  = sum((bsxfun(@times,derShareDeltij,qweightrprods)),4);

% Derivative wrt the random coefficient on segment

part1               = multixvuSEG ./(1-Sigseg);
part2A              = sum((bsxfun(@times,numer1,part1)),1);
part2               = Sigseg * (part2A ./denom1);
part2(isnan(part2)) = 0;

part3A                  = (1-Sigseg) * ((denom1 .^ (-Sigseg)) .* (part2A));
part3A(isnan(part3A))   = 0;
part3B                  = sum(part3A,2);
part3                   = part3B ./ denom2;

derShareRC1    = bsxfun(@minus,part1,part2);
derShareRC2    = sij .* bsxfun(@minus,derShareRC1,part3);
derShareRCSEG  = sum((bsxfun(@times,derShareRC2,qweightrprods)),4);

% Derivative wrt the random coefficient on continuous variable

part1               = multixvuRC ./(1-Sigseg);
part2A              = sum((bsxfun(@times,numer1,part1)),1);
part2               = Sigseg * (part2A ./denom1);
part2(isnan(part2)) = 0;

part3A                  = (1-Sigseg) * ((denom1 .^ (-Sigseg)) .* (part2A));
part3A(isnan(part3A))   = 0;
part3B                  = sum(part3A,2);
part3                   = part3B ./ denom2;

derShareRC1    = bsxfun(@minus,part1,part2);
derShareRC2    = sij .* bsxfun(@minus,derShareRC1,part3);
derShareRC     = sum((bsxfun(@times,derShareRC2,qweightrprods)),4);


% Derivative wrt the sigma segment

part1          = (bsxfun(@plus,delta,(xvuSEG.*RCSeg + xvuRC.*RCCont))) ./ ((1-Sigseg)^2);
part1          = reshape(part1,[prods,1,nmkt,Nnodes]);
part1          = bsxfun(@times,part1,multidumseg);

part2A                  = sum((part1 .* numer1),1);
part2B                  = bsxfun(@rdivide,part2A,denom1);
part2B(isnan(part2B))   = 0;
part2                   = bsxfun(@minus,(- log(denom1)),(Sigseg *part2B));
part2(isinf(part2))     = 0;

part3A                  = bsxfun(@plus,(- log(denom1)),((1-Sigseg) *part2B));
part3B                  = sum(((denom1.^(1-Sigseg)) .* part3A),2);
part3B(isnan(part3B))   = 0;
part3                   = part3B ./ denom2;
part3(isnan(part3))     = 0;

derShareSigseg1 = bsxfun(@minus,part2,part3);
derShareSigseg2 = sij .* bsxfun(@plus,part1,derShareSigseg1);
derShareSigseg  = sum((bsxfun(@times,derShareSigseg2,qweightrprods)),4);

if numel(thetaNL) == 2
    derShareTheta   = [sum(derShareRC,2),sum(derShareSigseg,2)];
    derDeltaTheta   = multiprod((-multinv(derShareDelta)),derShareTheta);
    derDeltaTheta   = [ reshape(derDeltaTheta(:,1,:,:),[Data.nobs 1]) ,...
                    reshape(derDeltaTheta(:,2,:,:),[Data.nobs 1])]; % 2D market shares
else
    
derShareTheta   = [ sum(derShareRCSEG,2) , sum(derShareRC,2) , sum(derShareSigseg,2)];
derDeltaTheta   = multiprod((-multinv(derShareDelta)),derShareTheta);
derDeltaTheta   = [ reshape(derDeltaTheta(:,1,:,:),[Data.nobs 1]) ,...
                    reshape(derDeltaTheta(:,2,:,:),[Data.nobs 1]),...
                    reshape(derDeltaTheta(:,3,:,:),[Data.nobs 1])]; % 2D market shares
end
               
f = derDeltaTheta;